home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / c_check2d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  5.5 KB  |  249 lines
  1. #include <stdio.h>
  2. #include <sys/time.h>
  3.  
  4. #include "fft.h"
  5. #include "constant.h"
  6.  
  7. /*                        */
  8. /*    Precision dependant declarations    */
  9. /*                        */
  10. #ifdef    ZOMPLEX
  11.  
  12. #define TOLERANCE   DTOLERANCE
  13. typedef        double    this_half;
  14. typedef        zomplex    this_type;
  15.  
  16. #define        THIS_SUM    dsum_
  17. #define        THIS_GEN    zgen_
  18. #define        THIS_FT    zft2d_
  19. #define        THIS_SCAL    zscal2d
  20. #define        THIS_FFTI    zfft2di
  21. #define        THIS_FFT    zfft2d
  22.  
  23. #endif
  24. #ifdef    COMPLEX
  25.  
  26. typedef        float        this_half;
  27. typedef        complex        this_type;
  28.  
  29. #define TOLERANCE   STOLERANCE
  30.  
  31. #define        THIS_SUM    ssum_
  32. #define        THIS_GEN    cgen_
  33. #define        THIS_FT    cft2d_
  34. #define        THIS_SCAL    cscal2d
  35. #define        THIS_FFTI    cfft2di
  36. #define        THIS_FFT    cfft2d
  37.  
  38. #endif
  39.  
  40. /*                        */
  41.  
  42. this_half THIS_SUM();
  43.  
  44. void inimat_();
  45. void GetArguments();
  46. void get_values();
  47.  
  48. int size_x, ldx, size_y, n_trials;
  49. this_type *pa, *pb, *pRef, *pSave;
  50.  
  51. int speedup = 0;
  52.  
  53. main(argc,argv)
  54. int argc;
  55. char *argv[];
  56. {
  57.     int i, j, k, n_total, is_wrong, nerr;
  58.     double x, y, dx, dy, emax, sum, err, maxerr;
  59.     this_half    t;
  60.     this_type *ptr;
  61.  
  62. /* ******************************************************* */
  63.     GetArguments( argc, argv);
  64. /* ******************************************************* */
  65.  
  66.     srandom( (123*getpid()) | 0x01);
  67.  
  68.     for( ; n_trials > 0; n_trials --) { 
  69.     get_values();
  70.  
  71.     printf("\n");
  72.     printf( "%d <%3d,%3d>", n_trials, size_x, size_y);
  73.     fflush(stdout);
  74.  
  75.     n_total = ((size_x+1)*(size_y+1));
  76.     ldx = size_x + 1;
  77.     pa = (this_type *)malloc
  78.         ( (3 * n_total + size_x + FACTOR_SPACE + size_y + FACTOR_SPACE) * sizeof(this_type));
  79.     if( !(pa) ) {
  80.         fprintf( stderr, "Could not allocate ... Exiting\n");
  81.         exit (-1);
  82.     }
  83.     pb = pa + n_total;
  84.     pRef = pb + n_total;
  85.     pSave = pRef + n_total;
  86.  
  87. /* *******************************************************
  88.     Initialize arrays
  89. ******************************************************* */
  90.     THIS_GEN(pRef, &n_total);
  91.     bcopy( pRef, pa, n_total*sizeof(this_type));
  92.     bcopy( pRef, pb, n_total*sizeof(this_type));
  93.  
  94. /* *******************************************************
  95.     direct Fourier Transform
  96. ******************************************************* */
  97.         j = -1;
  98.     THIS_FT( &j, &size_x, &size_y, pb, &ldx);
  99.     pSave = THIS_FFTI( size_x, size_y, pSave);
  100.     bcopy( pRef, pa, n_total*sizeof(this_type));
  101.     is_wrong = THIS_FFT( -1, size_x, size_y, pa, ldx, pSave);
  102.  
  103.     emax = TOLERANCE*(size_x * size_y);
  104.     for(i = 0, nerr=0 ; i < n_total ; i ++) {
  105.         x = pa[i].re - pb[i].re;
  106.         y = pa[i].im - pb[i].im;
  107.         if( (t= x*x + y*y) > (emax)) {
  108.         nerr++;
  109.         }
  110.     }
  111.     if( !(nerr)){
  112.         fprintf( stderr, ".");
  113.     }
  114.     else { 
  115.         fprintf( stderr, "@%d@", nerr);
  116.     }
  117.     i = 2;
  118.     x = y = 0;
  119.     for( j= 0, ptr = pRef ; j < size_y ; j++, ptr += ldx) { 
  120.         x +=  THIS_SUM( &size_x, ptr, &i);
  121.         y +=  THIS_SUM( &size_x, ((char *)ptr) + sizeof(this_half), &i);
  122.     }
  123.     dx = x - pa[0].re;
  124.     dy = y - pa[0].im;
  125.     if( (t= dx *dx + dy*dy) > TOLERANCE){
  126.         fprintf( stderr, 
  127.         "\nError direct FFT(%d,%d) at index #0 : (%f,%f)<->(%f,%f)error(%f,%f)",
  128.         size_x, size_y, pa[0].re, pa[0].im,x,y, dx, dy);
  129.     }
  130.     dx = x - pb[0].re;
  131.     dy = y - pb[0].im;
  132.     if( (t= dx *dx + dy*dy) > TOLERANCE){
  133.         fprintf( stderr, 
  134.         "\nError direct DFT(%d,%d) at index #0 : (%f,%f)<->(%f,%f)error(%f,%f)",
  135.         size_x, size_y, pb[0].re, pb[0].im,x,y, dx, dy);
  136.     }
  137. /* *******************************************************
  138.     inverse Fourier Transform
  139. ******************************************************* */
  140.     is_wrong = THIS_FFT( 1, size_x, size_y, pa, ldx, pSave);
  141.  
  142.     i = 2;
  143.     x = y = 0;
  144.     for( j= 0, ptr = pb ; j < size_y ; j++, ptr += ldx) { 
  145.         x +=  THIS_SUM( &size_x, ptr, &i);
  146.         y +=  THIS_SUM( &size_x, ((char *)ptr) + sizeof(this_half), &i);
  147.     }
  148.     dx = x - pa[0].re;
  149.     dy = y - pa[0].im;
  150.     nerr = 0;
  151.     emax = size_x * size_y * TOLERANCE;
  152.     if( (t= dx *dx + dy*dy) > TOLERANCE){
  153.         fprintf( stderr, "#");
  154.         nerr ++;
  155.     }
  156.     if( !(nerr)){
  157.         fprintf( stderr, ".");
  158.     }
  159.  
  160.     t = 1. / (double)(size_x * size_y);
  161.  
  162.     THIS_SCAL( size_x, size_y, t, pa, ldx);
  163.  
  164.     emax = TOLERANCE;
  165.     emax = emax * emax;
  166.  
  167.     sum = 0.;
  168.     err= 0.;
  169.     nerr= 0;
  170.     maxerr= 0.;
  171.     for(i = 0 ; i < n_total ; i ++) {
  172.         sum += (pRef[i].re * pRef[i].re) + (pRef[i].im * pRef[i].im);
  173.         x = pa[i].re - pRef[i].re;
  174.         y = pa[i].im - pRef[i].im;
  175.         if( (t= x*x + y*y) > (emax))
  176.         nerr++;
  177.         err += t;
  178.         if( t > maxerr) maxerr = t;
  179.     }
  180.     err = sqrt(err / (double)(size_x*size_y));
  181.     sum = sqrt(sum / (double)(size_x*size_y));
  182.     maxerr = sqrt(maxerr);
  183.     printf(": avg:(%8.3g /%8.3g)= %8.3g <> maxerr: (%8.3g /%8.3g)= %8.3g", 
  184.     err, sum, err/sum, maxerr, sum, maxerr/sum);
  185.     if(nerr){
  186.         printf("\n%d errors detected \n", nerr);
  187.         exit(-2);
  188.     }
  189.     free(pa);
  190.     }
  191.     printf("\n");
  192.     return(0);
  193. }
  194.  
  195. int is_random;
  196.  
  197. void GetArguments( argc, argv)
  198. int argc;
  199. char *argv[];
  200. {
  201.     int i, j, k;
  202.     int nerror = 0;
  203.  
  204. #define ON    1
  205.  
  206.     n_trials = DEF_TRIALS;
  207.     is_random = 1;
  208.  
  209. /* ******************************************************* */
  210.     for ( i = 1 ; (i < argc) && (nerror != ON) ; i ++ ) {
  211.     if( argv[i][0] == '-') {
  212.         switch ( argv[i][1]) {
  213.         case 'i' :
  214.         case 'I' :  
  215.                 is_random = 0;
  216.                 break;
  217.         default  :  nerror = ON;
  218.         }
  219.     }
  220.     else {
  221.         if( i+1 > argc)
  222.         nerror = ON;
  223.         else { 
  224.         n_trials = atoi( argv[i]);
  225.         }
  226.     }
  227.     }
  228.     if( nerror == ON) {
  229.     fprintf( stderr, 
  230.         "Usage : %s [-p(arallel)] [n_trials]\n", argv[0]);
  231.     exit(-1);
  232.     }
  233. }
  234.  
  235. void get_values()
  236. {
  237.   if( is_random){ 
  238.     size_x = random() % MAX_2D + 1;
  239.     size_y = random() % MAX_2D + 1;
  240.     if( !(n_trials % 5))
  241.     size_y = size_x;
  242.   } else{
  243.     printf( "Enter SIZE_X : ");
  244.     scanf("%d", &size_x);
  245.     printf( "Enter SIZE_Y : ");
  246.     scanf("%d", &size_y);
  247.   }
  248. }
  249.